
# ------------------------------------------------------------------------------------------
# # Julia很快
#
# 基准测试程序数值（benchmarks）通常被用来比较编程语言的性能。这些benchmarks可能引起很长的讨论，首先是作为基准测试的到底是什么，再有就是究竟是什么导致不同语言间的
# 差异。这些简单的问题有时比你乍想之下复杂得多。
#
# 这个notebook目的在于让你实现一个简单的基准测试。通过这个notebook可以看到作者的配有四核Macbook Pro表现如何，或者你可以自己跑跑看。
#
# (材料来自于MIT的Steven Johnson的一场演讲:
# https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-
# registers.ipynb.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # 大纲
#
# - 定义sum函数
# - 实现并测试sum函数在...
#     - C (自定义)
#     - C (使用-ffast-math的自定义)
#     - python (内置)
#     - python (numpy)
#     - python (自定义)
#     - Julia (内置)
#     - Julia (自定义)
#     - Julia (使用SIMD的自定义)
# - benchmarks总结
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # `sum`: 一个很容易理解的函数
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 考虑**求和**函数`sum(a)`，它是用来计算
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# 其中$n$是`a`的长度。
# ------------------------------------------------------------------------------------------

a = rand(10^7)  # 在[0,1)上均匀分布的随机数组成的1维向量

sum(a)

# ------------------------------------------------------------------------------------------
# 预期结果应该在0.5 * 10^7左右，因为所有元素的平均值应该是0.5左右
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # 几种语言下的几种基准测试方法
# ------------------------------------------------------------------------------------------

@time sum(a)

@time sum(a)

@time sum(a)

# ------------------------------------------------------------------------------------------
# `@time`宏会产生受干扰的结果，所以不是做基准测试的最佳选择！
#
# 好在Julia有个`BenchmarkTools.jl`包来简单并准确地进行基准测试：
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add("BenchmarkTools")

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# #  1. C语言
#
# C语言对人类不友好但是对机器友好，经常作为一个黄金标准。在C语言的两倍以内通常来说就比较理想了。然而，即使仅考虑C语言，不同的C语言程序员也可能采用或不采用不同的优化方式。
#
# 原作者不会C语言，所以他看不懂下面这个cell中的代码，但是他知道Julia回话中可以编译、运行C语言代码。注意`"""`所包围的是个多行字符串。
# ------------------------------------------------------------------------------------------

using Libdl
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()  # 创建一个临时文件


# 通过将C_code导入gcc来编译成共享链接库
# (只有装了gcc才行):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# 定义一个调用这个C语言函数的Julia函数：
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a)  # 输入\approx后按<TAB>键得到≈符号

c_sum(a) - sum(a)  

≈  # `isapprox` 函数的别名

?isapprox

# ------------------------------------------------------------------------------------------
# 现在我们可以在Julia中直接对C语言代码进行基准测试了：
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: 最快$(minimum(c_bench.times) / 1e6) 毫秒")

d = Dict()  # 字典，也就是关联数组（键值对）
d["C"] = minimum(c_bench.times) / 1e6  # 毫秒
d

using Plots
gr()

using Statistics  # 导入统计包以计算标准差
t = c_bench.times / 1e6  # 毫秒
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# # 2. 使用-ffast-math的C语言
#
# 如果我们允许C重新安排浮点操作，它就会使用SIMD指令(single instruction, multiple data)进行矢量化。
# ------------------------------------------------------------------------------------------

const Clib_fastmath = tempname()   # 创建临时文件

# 同上但是增加了-ffast-math标志
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# 定义一个调用这个C语言函数的Julia函数：
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_fastmath_bench = @benchmark $c_sum_fastmath($a)

d["C -ffast-math"] = minimum(c_fastmath_bench.times) / 1e6  # 毫秒

# ------------------------------------------------------------------------------------------
# # 3. Python的内置`sum`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `PyCall`包给Julia提供了Python接口：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("PyCall")
using PyCall

# 获得Python的内置"sum"函数：
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($a)

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 4. Python: `numpy`
#
# ## 条件允许的情况下利用硬件"SIMD"
#
# `numpy`是一个Python的优化过的C语言库。
# 在Julia版Python中的安装方法为：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("Conda")
using Conda

Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]

py_numpy_bench = @benchmark $numpy_sum($a)

numpy_sum(a)

numpy_sum(a) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 5. Python（自定义）
# ------------------------------------------------------------------------------------------

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($a)

sum_py(a)

sum_py(a) ≈ sum(a)

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 6. Julia (内置)
#
# ## 是直接用Julia写的，不是C！
# ------------------------------------------------------------------------------------------

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 7. Julia (自定义)
# ------------------------------------------------------------------------------------------

function mysum(A)   
    s = 0.0  # s = zero(eltype(a))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 8. Julia (使用simd的自定义)
# ------------------------------------------------------------------------------------------

function mysum_simd(A)   
    s = 0.0  # s = zero(eltype(A))
    @simd for a in A
        s += a
    end
    s
end

j_bench_hand_simd = @benchmark mysum_simd($a)

mysum_simd(a)

d["Julia hand-written simd"] = minimum(j_bench_hand_simd.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 总结
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end
